package scu.maqiang.elastodynamics;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class ElastoDynamics1DL21 {
	public static void main(String[] args) {
		Mesh1L2 mesh = new Mesh1L2();
		mesh.line(0, 1, 128);
		
		FES1L21 fs = new FES1L21(mesh);
		SRMatrix K = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] u = new double[fs.GetNdof()];
		double[] v = u.clone();
		double[] a = u.clone();
		double[] v0 = u.clone();
		double[] a0 = u.clone();
		
		double time = 0.0, dt = 0.005, endTime = 10.0;
		double alpha = 0.25, delta = 0.5;
		double c0 = 1.0 / (alpha * dt * dt);
		double c1 = delta / (alpha * dt);
		double c2 = 1.00 / (alpha * dt);
		double c3 = 0.50 / alpha - 1.0;
		double c4 = delta / alpha - 1.0;
		double c5 = 0.5 * dt * (delta / alpha - 2.0);
		double c6 = dt * (1 - delta);
		double c7 = delta * dt;
		
		fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, K);
		fs.assembleMass(new double[]{1.0}, BVPType.LUMP_MASS, M);
		System.out.println(M);
		SRMatrix L = K.add(c0, M);
		double[] uva = u.clone();
		//
		ScalarFunc initFunc =  (x, label, param) -> Math.sin(0.5 * Math.PI * x[0]); 
		double[] u0 = fs.valueFromFunc(initFunc, null);
		System.out.println(MVO.toString(u0));
		double[] Ku0 = K.mul(u0);
		K.sort();
		System.out.println(K);
		NewIterSSolver mSolver = new NewIterSSolver(M);
		for(int i = 0; i <fs.GetNdof(); i++) {
			RHS[i] = RHS[i] - Ku0[i];
		}
		mSolver.PCGSSOR(RHS, a0, 1.5, 1);
		Arrays.setAll(a, i -> a0[i]);
		
		for(int i = 0; i < fs.GetNdof(); i++) {
			uva[i] = c0 * u0[i] + c2 * v0[i] + c3 * a0[i];
		}
		M.mul(uva, RHS);
		String outputFile = "FreeVibration1D.dat";
		time = 0.0;
		mesh.toTecplot(outputFile, time, u0);
		fs.applyBC_MBN(L, Direct.All, 1);
//		L.sort();
//		//System.out.println(L);
//		System.out.println(L.getNNZ());
		
		NewIterSSolver solver = new NewIterSSolver(L);
//		ScalarCoefFunc BCValueFunc = (xy, llabel, param) -> 0.1 * Math.sin(8 * Math.PI * param[0]);
		
		//myParam[0] = 0.005;
		//fs.applyBC_MBN(RHS, Direct.X, new RegionFunc[] {BCValueFunc}, myParam, exteriorLabel);
		//for(int i = 0; i < fs.GetNdof(); i++) {
		//	System.out.println(RHS[i]);
		//}
//		ScalarCoefFunc[] bdRegionFunc = new ScalarCoefFunc[] {BCValueFunc};

		double[] myParam = new double[1];
		while(time <= endTime) {
			time += dt;
			myParam[0] = time;
			System.out.println("time: " + time);
			fs.applyBC_MBN(RHS, Direct.All, 0.0, 1);
			//System.out.println(Arrays.toString(RHS));
			solver.PCGSSOR(RHS, u, 1.5, 1);
			Arrays.setAll(a, i -> c0 * (u[i] - u0[i]) - c2 * v0[i] - c3 * a0[i]);
			Arrays.setAll(v, i -> v0[i] + c6 * a0[i] + c7 * a[i]);
			Arrays.setAll(uva, i -> c0 * u[i] + c2 * v[i] + c3 * a[i]);
			M.mul(uva, RHS);
			Arrays.setAll(u0, i -> u[i]);
			Arrays.setAll(v0, i -> v[i]);
			Arrays.setAll(a0, i -> a[i]);
			mesh.toTecplot(outputFile, time, u0);
		}		
	}
}
